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1 Introduction 

On the nanometer length scale, continuum approaches like the Navier-Stokes equa- 
tion break down fl]. Therefore, the study of nanoscopic transport processes re- 
quires a molecular point of view and preferably the application of molecular dy- 
namics (MD) simulation. In the past, MD could be applied to small systems with 
a few thousand particles only, due to the low capacity of computing equipment. 
Consequently, a large gap between MD simulation results on the one hand and ex- 
perimental results as well as calculations based on continuum methods was present. 

The constant increase in available computational power is eliminating this bar- 
rier, and the characteristic length of the systems accessible to MD simulation ap- 
proaches micrometers. However, this can only be realized by laying an emphasis 
on the simplicity of the molecular models and an efficient implementation, suitable 
for massively parallel processing. 

Due to their anisotropy, nanostructures containing graphite and carbon nan- 
otubes are of particular interest. The present work deals with the flow behavior of 
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liquid methane, modeled by the truncated and shifted Lennard-Jones (LJ/TS) po- 
tential (2), confined between graphite walls. It covers nanoscopic Poiseuille flow 
up to a channel width of 0.135 /jm. 

2 Simulation method and scalability 

With the potential parameters a f = 3.7241 A and e^/k-g = 175.06 K, the LJ/TS 
potential is known to be an accurate model for methane that covers the thermody- 
namic properties of the fluid quantitatively Q. Following a widespread approach 
introduced by Battezzati et al. [4], the interaction between methane molecules and 
the carbon atoms begin part of a graphite surface can also be be described by a 
Lennard-Jones potential. The present study applies the LJ/TS potential to the in- 
teraction between methane molecules and carbon atoms as well, using the poten- 
tial parameters proposed for graphite by Wang et al. |5], <r w = 3.3264 A and e w 
= 0.00188 eV. The unlike interaction parameters <7f w and £f w , acting between 
methane and carbon, were determined according to the Lorentz-Berthelot mixing 
rule. 

Carbon and silicon structures as well as ceramics can well be represented by the 
Tersoff potential J6JI3, which permits to predict many properties of these systems 
with a good accuracy (H. For graphite, however, the bond length corresponding 
to the Tersoff potential (1.461 A), cf. Kelires [9], deviates considerably from the 
actual value of 1.421 A iflOl . Therefore, the relevant Tersoff potential parameters 
for the wall model were rescaled in the present study. 

The present version of the employed LS1/Mardyn MD simulator uses a spa- 
tial domain decomposition with equally sized cuboid subdomains and a cartesian 
topology based on linked cells [11 J. Often the best solution is an 'isotropic' de- 
composition that minimizes the surface to volume ratio of the spatial subdomains. 
For the simulation of homogeneous systems, this approach is quite efficient f!2l . 
That is underlined by the weak and strong scaling behavior of LS1/Mardyn for 
typical configurations, shown in Fig.[T](left), in cases where supercritical methane 
('fluid') at p = 10 mol/1 and solid graphite ('wall') were considered with a system 
size of up to 4,800,000 interaction sites, representing the same number of carbon 
atoms and methane molecules here. Graphite simulations, where only the carbon 
wall atoms are regarded, scale particularly well, due to a favorable relation of the 
delay produced by communication between processes to the concurrent parts, i.e., 
the actual intermolecular interaction computation, which is much more expensive 
for the Tersoff potential than the LJ/TS potential. 

The simulation of the regarded combined systems, containing both fluid and 
solid interaction sites, is better handled by a channel geometry based decomposi- 
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Figure 1: Left: total CPU time, i.e., execution time multiplied with the number of 
parallel processes, per time step and interaction site for weak scaling with 3,000 
(dashed lines / circles) and 32,000 (dotted lines / triangles) interaction sites per pro- 
cess as well as strong scaling with 450,000 interaction sites (solid lines / bullets), 
using isotropic spatial domain decomposition. Right: speedup, i.e., sequential ex- 
ecution time divided by parallel execution time, for a system of liquid methane 
between graphite walls with 650,000 interaction sites, where isotropic (circles) 
and channel geometry based (triangles) spatial domain decomposition was used; 
the solid line represents optimal speedup. 

tion scheme, where an approximately equal portion of the wall and a part of the 
fluid is assigned to each process, cf. Fig. Q] (right). In the general case, where spa- 
tial non-uniformities do not match any cartesian grid, a flexible topology has to be 
used. An approach based on /c-dimensional trees |[T3l[T4l . implemented in a ver- 
sion of LS1/Mardyn, showed clearly improved results with respect to the scaling 
of inhomogeneous systems. 

3 MD flow simulation of liquid methane 

The LS1/Mardyn MD program was used to conduct a series of Poiseuille flow 
simulations with liquid methane in the canonical ensemble. The flow was induced 
by an external gravitation-like acceleration acting on all fluid molecules. Anal- 
ogously, Couette flow was simulated by accelerating only the wall atoms. A PI 
controller 

r 2 a z (t)=v z -2v z (t)+v z {t-r'), (1) 
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Figure 2: Left: velocity profile (top) and density profile (bottom) for Poiseuille 
(solid lines) and Couette (dotted lines) flow of liquid methane at T = 166.3 K with 
a channel width of h = 8 nm and a characteristic flow velocity of v z = 50 m/s. Right: 
pressure drop — Ap in terms of v z and the channel length Az (top) as well as slip 
velocity in terms of v z (bottom), for Poiseuille flow of saturated liquid methane at 
a temperature of T = 166.3 K and average velocities v z of 10 m/s (circles) and 30 
m/s (bullets), in dependence of the channel width; solid line: Darcy's law. 



was applied to the acceleration a z for a flow in z direction with a characteristic 
velocity of v z , where v z (t) is the velocity at a given time, r was on the order of 1 
- 100 ps and r' on the order of 0. 1 - 10 ps. Velocity and density profiles, cf. Fig. |2] 
(left), and the average acceleration were extracted from the simulations. 

It is known for Poiseuille flow with channel widths h below 2 nm that the slip 
velocity can reach values above 0.99 v z , a regime which was studied by Sokhan et 
al. lfl5l for methane in carbon nanotubes. Velocity profiles with a strong influence 
of boundary slip were also found in some of the present simulations, cf. Fig. |2] 
(left). For channel widths between h = 20 and 50 nm, the boundary slip undergoes 
a transition, as shown in Fig. |2] (right). For h down to molecular length scales, the 
pressure drop — Ap is approximately proportional to the average velocity v z and 
inversely proportional to cross-sectionional area of the channel, corresponding to 
cf. Fig. |2] (right), in agreement with Darcy's law. 
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4 Conclusion 



The scalability of the LS1/Mardyn program, optimized for massively parallel 
MD simulations of simple fluids interacting with carbon nanostructures, was as- 
sessed and found to be satisfactory. MD simulations of methane confined be- 
tween graphite walls with up to 4,800,000 interaction sites, i.e., carbon atoms and 
methane molecules, were conducted to demonstrate the viability of the program. 

The channel width was varied to include both the boundary-dominated regime 
and the transition to the continuum regime. This proves that MD can be used today 
to cover the entire range of characteristic lengths for which continuum methods 
fail. The simulation results show that the transition between both regimes occurs 
in a relatively narrow region, between h = 20 and 50 nm in the present case. For a 
flow velocity up to 30 m/s, it was confirmed for methane in graphite channels that 
Darcy's law applies on both sides of this transition. 
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